{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import numpy.linalg as la\n",
    "\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.preprocessing import scale\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import cm as cm\n",
    "import seaborn as sns\n",
    "sns.set(font_scale=2)\n",
    "plt.style.use('seaborn-whitegrid')\n",
    "%matplotlib inline\n",
    "\n",
    "from io import StringIO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Principal Component Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In machine learning and data science, Principal Component Analysis (PCA) is a method of analysing datasets to obtain an orthogonal basis that best captures most of the variance of our data.  It is often used to remove extra features or dimensions in large-dimensional datasets, because in a similar vein to SVD, PCA will give us a set of axes with both large contributions and small contributions to our data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 1: Breast cancer dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Importing the data set for characteristics of tumor cells (this is the same dataset used for the MP)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = [\"radius\", \"texture\", \"perimeter\", \"area\",\n",
    "          \"smoothness\", \"compactness\", \"concavity\",\n",
    "          \"concave points\", \"symmetry\", \"fractal dimension\"];\n",
    "stats = [\"(mean)\", \"(stderr)\", \"(worst)\"]\n",
    "labels = [\"patient ID\", \"Malignant/Benign\"]\n",
    "\n",
    "for p in params:\n",
    "    for s in stats:\n",
    "        labels.append(p + \" \" + s)\n",
    "\n",
    "tumor_data = pd.io.parsers.read_csv(\"additional_files/breast-cancer-train.dat\",header=None,names=labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tumor_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1) Creating a smaller subset of points:\n",
    "Selecting a subset of the data for better visualization and understanding of the method. We will start with six patients and only two of the features:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_data = pd.DataFrame(tumor_data[[\"Malignant/Benign\", 'smoothness (mean)', 'radius (mean)']][272:278])\n",
    "new_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g1 = sns.lmplot('smoothness (mean)', 'radius (mean)', new_data,  hue=\"Malignant/Benign\", scatter_kws={\"s\": 180}, fit_reg=False, height=8)\n",
    "ax = g1.axes[0,0]\n",
    "ax.axis('equal')\n",
    "\n",
    "# This code snippet is plotting the labels\n",
    "for i in range(272,278):\n",
    "    x = new_data['smoothness (mean)'][i] + 0.1\n",
    "    y = new_data['radius (mean)'][i] + 0.1\n",
    "    ax.text(x,y,str(i),horizontalalignment='left',size='medium', color='black', weight='semibold', fontsize=16)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2) Shift the dataset to center the data:\n",
    "The first step is to determine the \"center\" of the dataset (the mean value of each feature):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_smooth = new_data['smoothness (mean)'].mean()\n",
    "mean_radius = new_data['radius (mean)'].mean()\n",
    "\n",
    "print(mean_smooth,mean_radius)\n",
    "\n",
    "print(new_data['smoothness (mean)'].std(),new_data['radius (mean)'].std())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g1 = sns.lmplot('smoothness (mean)', 'radius (mean)', new_data,  hue=\"Malignant/Benign\", scatter_kws={\"s\": 180}, fit_reg=False, height=8)\n",
    "ax = g1.axes[0,0]\n",
    "plt.xlim(-2,24)\n",
    "plt.ylim(-2,24)\n",
    "for i in range(272,278):\n",
    "    x = new_data['smoothness (mean)'][i] + 0.1\n",
    "    y = new_data['radius (mean)'][i] + 0.1\n",
    "    ax.text(x,y,str(i),horizontalalignment='left',size='medium', color='black', weight='semibold', fontsize=16)\n",
    "        \n",
    "ax.scatter(mean_smooth,mean_radius,  s=180, c='r', marker=(5, 2))\n",
    "ax.axhline(y=0, color='k')\n",
    "ax.axvline(x=0, color='k')\n",
    "ax.axhline(y=mean_radius, color='m',linestyle='--')\n",
    "ax.axvline(x=mean_smooth, color='m',linestyle='--')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And then we \"center\" the dataset, such that each feature has zero mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_data['smoothness (zero mean)'] = new_data['smoothness (mean)'] - new_data['smoothness (mean)'].mean()\n",
    "new_data['radius (zero mean)'] = new_data['radius (mean)'] - new_data['radius (mean)'].mean()\n",
    "\n",
    "print(new_data['smoothness (zero mean)'].mean())\n",
    "print(new_data['radius (zero mean)'].mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_data['smoothness (zero mean)']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the centered data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g1 = sns.lmplot('smoothness (zero mean)', 'radius (zero mean)', new_data, hue=\"Malignant/Benign\", scatter_kws={\"s\": 180}, fit_reg=False, height=8)\n",
    "ax = g1.axes[0,0]\n",
    "ax.axis('equal')\n",
    "\n",
    "for i in range(272,278):\n",
    "    x = new_data['smoothness (zero mean)'][i] + 0.1\n",
    "    y = new_data['radius (zero mean)'][i] + 0.1\n",
    "    ax.text(x,y,str(i),horizontalalignment='left',size='medium', color='black', weight='semibold', fontsize=16)\n",
    "\n",
    "ax.scatter(0,0,  s=200, c='r', marker=(5, 2))\n",
    "\n",
    "ax.axhline(y=0, color='k')\n",
    "ax.axvline(x=0, color='k')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3) Get covariance matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PCA wants to find the directions of maximum variance. For that, we will need to first define the covariance matrix:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$ M $: total number of data points\n",
    "\n",
    "$ N $ : number of features\n",
    "\n",
    "$Cov({\\bf A}) = \\frac{1}{M-1} {\\bf A}^T {\\bf A} $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# centered data\n",
    "A = new_data[['smoothness (zero mean)', 'radius (zero mean)']]\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M,N = A.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using the definition\n",
    "cov_matrix = (1/(M-1))*A.T@A\n",
    "print(cov_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Or using python built-in function\n",
    "A.cov()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The variances are on the diagonal (co-variance of a variable with itself), and the sum of the 2 values is the overall variability:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.diag(cov_matrix).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "7.53951/31.35945"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PCA replaces the original variables with new variables, called principal components, which are orthogonal (i.e. they have zero covariations) and have variances in decreasing order. To accomplish this, we will use the diagonalization of the covariance matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l,u = la.eig(cov_matrix)\n",
    "print(l)\n",
    "print(u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$cov({\\bf A}) = \\begin{bmatrix} -0.40237101 & -0.91547669 \\\\ -0.91547669 &  0.40237101 \\end{bmatrix} \\begin{bmatrix} 27.71798127 & 0 \\\\ 0 &  3.6414725  \\end{bmatrix} \\begin{bmatrix} -0.40237101 & -0.91547669 \\\\ -0.91547669 &  0.40237101 \\end{bmatrix}^T$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "27.717981/31.3594537"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the diagonal sum is still 31.359, which says that the two components account for all the variability.\n",
    "\n",
    "\n",
    "PCA finds, in the data space, the dimension (direction) with the largest variance out of the overall variance.\n",
    "\n",
    "\n",
    "In this example, if we reduce the dimension space to include only one variable, the first principal component 27.718, accounts for 88% of the variability"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hence, the largest eigenvalue of the covariance matrix corresponds to the largest variance of the dataset, and the associated eigenvector is the direction of maximum variance. For our example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g1 = sns.lmplot('smoothness (zero mean)', 'radius (zero mean)', new_data, hue=\"Malignant/Benign\", scatter_kws={\"s\": 180}, fit_reg=False, height=8)\n",
    "ax = g1.axes[0,0]\n",
    "ax.axis('equal')\n",
    "\n",
    "for i in range(272,278):\n",
    "    x = new_data['smoothness (zero mean)'][i] + 0.1\n",
    "    y = new_data['radius (zero mean)'][i] + 0.1\n",
    "    ax.text(x,y,str(i),horizontalalignment='left',size='medium', color='black', weight='semibold', fontsize=16)\n",
    "\n",
    "ax.scatter(0,0,  s=200, c='r', marker=(5, 2))\n",
    "\n",
    "ax.axhline(y=0, color='k')\n",
    "ax.axvline(x=0, color='k')\n",
    "\n",
    "s = 3\n",
    "\n",
    "J = 0 # choice of principal direction\n",
    "x = u[0,J]\n",
    "y = u[1,J]\n",
    "ax.arrow(0,0,s*x,s*y,color='black',head_width=0.1, head_length=0.1, fc='r', ec='r', lw=5)\n",
    "\n",
    "J = 1 # choice of principal direction\n",
    "x = u[0,J]\n",
    "y = u[1,J]\n",
    "ax.arrow(0,0,s*x,s*y,color='black',head_width=0.1, head_length=0.1, fc='m', ec='m', lw=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4) Singular value decomposition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We know that the eigenvectors of ${\\bf A}^T{\\bf A}$ are the right singular vectors of ${\\bf A}$, or the columns of ${\\bf V}$ from the SVD decomposition of ${\\bf A}$ (or the rows of V transpose). \n",
    "\n",
    "Hence, instead of having to calculate the covariance matrix and solve an eigenvalue problem, we will instead get the reduced form of the SVD!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "U, S, Vt = np.linalg.svd(A, full_matrices=False)\n",
    "\n",
    "# variances = eig(covariance) = singular values squared\n",
    "variances = S**2\n",
    "\n",
    "print(variances)\n",
    "\n",
    "# principal directions\n",
    "pc1_vec = Vt[0,:]\n",
    "pc2_vec = Vt[1,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we don't need to worry about the constant value from the covariance matrix that we are disregarding. The variance values change, but their proportionality remains:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "variances[0]/variances.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5) Plotting the principal directions using the singular right vectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g1 = sns.lmplot('smoothness (zero mean)', 'radius (zero mean)', new_data, hue=\"Malignant/Benign\", scatter_kws={\"s\": 180}, fit_reg=False, height=8)\n",
    "ax = g1.axes[0,0]\n",
    "\n",
    "for i in range(272,278):\n",
    "    x = new_data['smoothness (zero mean)'][i] + 0.1\n",
    "    y = new_data['radius (zero mean)'][i] + 0.1\n",
    "    ax.text(x,y,str(i),horizontalalignment='left',size='medium', color='black', weight='semibold', fontsize=16)\n",
    "\n",
    "ax.scatter(0,0,  s=200, c='r', marker=(5, 2))\n",
    "ax.axis('equal')\n",
    "\n",
    "ax.axhline(y=0, color='k')\n",
    "ax.axvline(x=0, color='k')\n",
    "\n",
    "s = 3\n",
    "\n",
    "x = pc1_vec[0]\n",
    "y = pc1_vec[1]\n",
    "ax.arrow(0,0,s*x,s*y,color='black',head_width=0.1, head_length=0.1, fc='r', ec='r', lw=5)\n",
    "\n",
    "x = pc2_vec[0]\n",
    "y = pc2_vec[1]\n",
    "ax.arrow(0,0,s*x,s*y,color='black',head_width=0.1, head_length=0.1, fc='m', ec='m', lw=5)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6) Cumulative explained variance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "variances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(pc1_vec)\n",
    "print(pc2_vec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, the largest variance is 138.6 and the direction of this principal component is given by the vector `pc1_vec`.\n",
    "\n",
    "The second largest variance is 18.2 and the direction of this principal component is given by the vector `pc2_vec`.\n",
    "\n",
    "In a general problem, we would have many principal components. How can we easily visualize these components and decide how many we will keep in our reduced feature space? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tot = sum(variances)\n",
    "var_exp = [(i / tot)*100 for i in variances]\n",
    "cum_var_exp = np.cumsum(var_exp)\n",
    "\n",
    "plt.bar(range(len(var_exp)),var_exp, align='center', label='individual explained variance')\n",
    "plt.step(range(len(var_exp)), cum_var_exp, 'r', where='mid', label='cumulative explained variance')\n",
    "plt.ylabel('Explained variance ratio')\n",
    "plt.xlabel('Principal components')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xstar=(A@Vt.T).values\n",
    "new_data['pc1'] = Xstar[:,0]  \n",
    "new_data['pc2'] = Xstar[:,1]  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g1 = sns.lmplot('pc1', 'pc2', new_data, hue=\"Malignant/Benign\", fit_reg=False, height=8, scatter_kws={\"s\": 180})\n",
    "ax = g1.axes[0,0]\n",
    "ax.axhline(y=0, color='k')\n",
    "ax.axvline(x=0, color='k')\n",
    "ax.axis('equal')\n",
    "\n",
    "\n",
    "for i in range(272,278):\n",
    "    x = new_data['pc1'][i] + 0.1\n",
    "    y = new_data['pc2'][i] + 0.1\n",
    "    ax.text(x,y,str(i),horizontalalignment='left',size='medium', color='black', weight='semibold', fontsize=16)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7) Complete dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But since we have only two features here, PCA is not really helping! Let's go back to original example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tumor_data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A_large = tumor_data.iloc[:,2:].values\n",
    "\n",
    "type(A_large)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A_large.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Center the mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A_large.std(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# X = (A_large - A_large.mean(axis=0))\n",
    "# print(X.mean(axis=0))\n",
    "\n",
    "X = (A_large - A_large.mean(axis=0))/A_large.std(axis=0)\n",
    "print(X.std(axis=0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u,s,vt = la.svd(X,full_matrices=False)\n",
    "\n",
    "variances = s**2\n",
    "\n",
    "tot = sum(variances)\n",
    "var_exp = [(i / tot)*100 for i in variances]\n",
    "cum_var_exp = np.cumsum(var_exp)\n",
    "\n",
    "plt.bar(range(len(var_exp)),var_exp, align='center', label='individual explained variance')\n",
    "plt.step(range(len(var_exp)), cum_var_exp, 'r', where='mid', label='cumulative explained variance')\n",
    "plt.ylabel('Explained variance ratio')\n",
    "plt.xlabel('Principal components')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cum_var_exp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Get the most important principal directions, and transform the original dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Vstar = vt[:3,:].T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xstar=(X@Vstar) # change of basis\n",
    "\n",
    "Xstar.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tumor_data_new = tumor_data.copy()\n",
    "\n",
    "tumor_data_new['pc1'] = Xstar[:,0]  \n",
    "tumor_data_new['pc2'] = Xstar[:,1]  \n",
    "tumor_data_new['pc3'] = Xstar[:,2]  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g1 = sns.lmplot('pc1', 'pc2', tumor_data_new, hue=\"Malignant/Benign\", fit_reg=False, height=8, scatter_kws={\"s\": 180})\n",
    "ax = g1.axes[0,0]\n",
    "ax.axhline(y=0, color='k')\n",
    "ax.axvline(x=0, color='k')\n",
    "ax.axis('equal')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mpl_toolkits.mplot3d import Axes3D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "threedee = plt.figure().gca(projection='3d')\n",
    "threedee.scatter(tumor_data_new['pc1'], tumor_data_new['pc2'], tumor_data_new['pc3'],c = tumor_data_new[\"Malignant/Benign\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot the weight for each feature in the first principal component:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "feature_names = tumor_data_new.columns[2:32]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "V = vt.T\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14,6))\n",
    "plt.bar(feature_names,V[:,0])\n",
    "plt.xticks(rotation=90);\n",
    "plt.title('importance of each attribute in ${\\\\bf p}_1$');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 2: Principal Components of FIFA Dataset\n",
    "We will be looking at the [FIFA 2018 Dataset](https://www.kaggle.com/thec03u5/fifa-18-demo-player-dataset/kernels).  While this is a video game, the developers strive to make their game as accurate as possible, so this data reflects the skills of the real-life players.\n",
    "\n",
    "Let's load the data frame using `pandas`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(\"additional_files/FIFA_2018.csv\",encoding = \"ISO-8859-1\",index_col = 0, low_memory = False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can take a brief look at the data by calling `df.head()`.  The first 34 columns are attributes that describe the behavior (e.g. aggression) or the skills (e.g. ball control), of each player.  The final columns show the player's position, name, nationality, and the club they play for.\n",
    "\n",
    "The four positions are forward (FWD), midfielder (MID), defender (DEF), and goalkeeper (GK)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A higher number signifies that an attribute is more prevalent for that player.  Looking at the above rankings, Player 0 (Christiano Ronaldo) has very good ball control and composure, but is not overly aggressive."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Correlation Matrix\n",
    "We can compute the correlation matrix for these variables across all players using a \"heatmap\".  Calling `df.corr()` provides this correlation matrix, and `seaborn.heatmap` will do the plotting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15,8))\n",
    "sns.heatmap(df.corr(),vmin=-1.0,vmax=1.0, linewidth=0.25, cmap='coolwarm');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This heatmap is dark red whenever two variables are positively correlated, and dark blue when they are negatively correlated.  For example, \"Sprint Speed\" and \"Acceleration\" are positively correlated.  \"Balance\" and \"Strength\" are negatively correlated however.\n",
    "\n",
    "Notice across the diagonal, all rectangles are dark red.  This is to be expected, as any variable is perfectly correlated with itself.\n",
    "\n",
    "Also notice that all Goal-Keeping skills are positively correlated with each other, but are negatively correlated with nearly all the other variables.  Maybe we can compress these into a single component/feature through principal component analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Principal Component analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that Principal Component Analysis (PCA) projects high-dimensional data into a low-dimensional representation by finding directions of maximal variance.\n",
    "\n",
    "Let's first create a new dataframe that includes only the attributes of each player (and not the last four columns of `df`). Store this new dataframe as a variable `X`. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = df.iloc[:,:-4].copy()\n",
    "X"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can get all the attribute names and store them as `labels` by using `.columns.values`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "labels = X.columns.values\n",
    "labels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To perform PCA, we first shift the data so that each attribute has zero mean, then compute the Singular Value Decomposition (SVD) of the resulting data matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create the data frame `A` where each attribute has zero mean.  Should we ensure each row has zero mean, or each column?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = X - X.mean()\n",
    "A.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now compute the SVD of the resulting matrix.  Make sure you compute the *reduced* SVD.\n",
    "\n",
    "Once you have computed the SVD, you can plot the fraction of explained variance for each singular value\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\sigma_i^2}{\\sum_{k=1}^r\\sigma_k^2} \\hspace{7mm} i = 1,2,\\dots,r\n",
    "\\end{equation}\n",
    "\n",
    "as well as the cumulative explained variance\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\sum_{k=1}^i \\sigma_k^2}{\\sum_{k=1}^r\\sigma_k^2} \\hspace{7mm} i = 1,2,\\dots,r\n",
    "\\end{equation}\n",
    "\n",
    "You can create a bar plot of the fraction of explained variance for each singular value using `plt.bar`, and a standard line plot for the cumulative explained variance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "U, S, Vt = np.linalg.svd(A, full_matrices = False)\n",
    "V = Vt.T\n",
    "\n",
    "variance = S**2\n",
    "sum_var = sum(variance)\n",
    "var_exp = [v/sum_var for v in variance]\n",
    "cum_var = np.cumsum(var_exp)\n",
    "\n",
    "plt.figure(figsize=(10,6))\n",
    "plt.bar(range(34),var_exp,label='individual explained variance')\n",
    "plt.plot(range(34),cum_var,'ro-', label='cumulative explained variance')\n",
    "plt.legend(loc=5)\n",
    "plt.xlabel(\"components\")\n",
    "plt.ylabel(\"% variance\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cum_var"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should see from the graph that the first principal component is responsible for nearly 60% of the variance, and the first two principal components have well over 70%.\n",
    "\n",
    "Recall from the SVD that $\\mathbf{A}\\mathbf{v}_i = \\sigma_i\\mathbf{u}_i$.  Writing the columns of $\\mathbf{A}$ as $\\mathbf{a}_k$, this means that:\n",
    "\n",
    "\\begin{equation}\n",
    "\\mathbf{v}_i^{(1)}\\begin{bmatrix}\\vdots \\\\ \\mathbf{a}_1 \\\\ \\vdots\\end{bmatrix} + \\mathbf{v}_i^{(2)}\\begin{bmatrix}\\vdots \\\\ \\mathbf{a}_2 \\\\ \\vdots\\end{bmatrix} + \\dots + \\mathbf{v}_i^{(n)}\\begin{bmatrix}\\vdots \\\\ \\mathbf{a}_n \\\\ \\vdots\\end{bmatrix} = \\sigma_i\\mathbf{u}_i\n",
    "\\end{equation}\n",
    "\n",
    "where $\\mathbf{v}_i^{(j)}$ is the $j$-th component of $\\mathbf{v}_i$. Thus if we define the principal components as \n",
    "\n",
    "$${\\bf p}_i = \\sigma_i\\mathbf{u}_i,$$ \n",
    "\n",
    "\n",
    "the $i$-th column of $\\mathbf{V}$ describes the projection of each attribute onto that principal direction."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can visualize the weight of each attribute to a given principal component by plotting the entries of the corresponding column of ${\\bf V}$. For example, the plot below illustrates the \"importance\" of each attribute to the first principal component (${\\bf p}_1$)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14,6))\n",
    "plt.bar(labels,V[:,0])\n",
    "plt.xticks(rotation=90);\n",
    "plt.title('importance of each attribute in ${\\\\bf p}_1$');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14,6))\n",
    "plt.bar(labels,V[:,1])\n",
    "plt.xticks(rotation=90);\n",
    "plt.title('importance of each attribute in ${\\\\bf p}_2$');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's add two new columns to the original dataframe `df`, with headers `pc1` and `pc2`. \n",
    "\n",
    "Use the expression above to evaluate the first two principal components ${\\bf p}_1$ and ${\\bf p}_2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# clear\n",
    "df['pc1'] = U[:,0]*S[0]\n",
    "df['pc2'] = U[:,1]*S[1]\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the data with these first two principal components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = sns.lmplot(x = \"pc1\", y = \"pc2\", data = df, hue = \"Position\", fit_reg=False, height=11, aspect=2, legend=True,\n",
    "           scatter_kws={'s':14,'alpha':0.5})\n",
    "ax = g.axes[0,0]\n",
    "ax.axvline(x=0,color='k', ls = '--')\n",
    "ax.axhline(y=0,color='k', ls = '--')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like the first principal axis determines whether a player is a goalkeeper or not.  We should double-check to make sure."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What are the attributes of $A$ that are most positively correlated with the first principal component? \n",
    "\n",
    "We can answer that by looking at the plot of coefficients above. Or we can do this in a systematic way, by sorting the entries of the column of ${\\bf V}$ and finding the ones with highest positive values. \n",
    "\n",
    "Find the first 5 attributes, and print their corresponding weights."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ind = np.argsort(V[:,0])\n",
    "print(ind)\n",
    "print(labels[ind[-5:]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see that all the goalkeeper attributes are positively correlated with the first principal component.  However, all other attributes, beginning with \"Strength\" are negatively correlated.  Try plotting the projection of \"GK reflexes\" onto the first two principal components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.where((labels == \"Finishing\") == True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.where((labels == \"GK reflexes\") == True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = sns.lmplot(x = \"pc1\", y = \"pc2\", data = df, hue = \"Position\", fit_reg=False, height=11, aspect=2, legend=True,\n",
    "            markers=[\"o\", \"x\",\"^\",\"s\"],palette=dict(FWD=\"g\", GK=\"orange\", MID=\"r\", DEF=\"m\"))\n",
    "ax = g.axes[0,0]\n",
    "ax.axvline(x=0,color='k', ls = '--')\n",
    "ax.axhline(y=0,color='k', ls = '--')\n",
    "\n",
    "s = 400   # this will scale the size of the arrow plot\n",
    "J = 15        # looking at the position \"GK reflexes\", corresponding to column 31\n",
    "x = V[J,0]    # projection of \"GK reflexes\" onto first principal component\n",
    "y = V[J,1]    # projection of \"GK reflexes\" onto second principal component\n",
    "\n",
    "# make an arrow from the origin to a point at (x,y)\n",
    "ax.arrow(0,0,s*x,s*y,color='black',width=1)\n",
    "ax.text(x*s*1.05,y*s*1.05,labels[J],fontsize=24)\n",
    "\n",
    "J = 9         # looking at the position \"GK reflexes\", corresponding to column 31\n",
    "x = V[J,0]    # projection of \"GK reflexes\" onto first principal component\n",
    "y = V[J,1]    # projection of \"GK reflexes\" onto second principal component\n",
    "\n",
    "# make an arrow from the origin to a point at (x,y)\n",
    "ax.arrow(0,0,s*x,s*y,color='black',width=1)\n",
    "ax.text(x*s*1.05,y*s*1.05,labels[J],fontsize=24)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you plot any other of the GK attributes, they will essentially overlap with GK reflexes.  Check that, by changing the variable `J` above to take the values (11,12,13,14)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make the same plot as above, but now take a look at other attributes. In the same figure, plot the projections for the attributes in columns [1,8,9,16,28,31]. \n",
    "\n",
    "Do you think the results make sense? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#clear\n",
    "\n",
    "g = sns.lmplot(x = \"pc1\", y = \"pc2\", data = df, hue = \"Position\", fit_reg=False, height=11, aspect=2, legend=True,\n",
    "            markers=[\"o\", \"x\",\"^\",\"s\"],palette=dict(FWD=\"g\", GK=\"orange\", MID=\"r\", DEF=\"m\"))\n",
    "ax = g.axes[0,0]\n",
    "ax.axvline(x=0,color='k', ls = '--')\n",
    "ax.axhline(y=0,color='k', ls = '--')\n",
    "\n",
    "s = 300   # this will scale the size of the arrow plot\n",
    "\n",
    "for J in [1,8,9,16,28,31]:\n",
    "\n",
    "    x = V[J,0]    \n",
    "    y = V[J,1]    \n",
    "\n",
    "    # make an arrow from the origin to a point at (x,y)\n",
    "    ax.arrow(0,0,s*x,s*y,color='black',width=1)\n",
    "    ax.text(x*s*1.5,y*s*1.1,labels[J],fontsize=24)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Remove data and re-do PCA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first principal component seems to mainly dictate whether a player is a goal-keeper or not.  To find out more about the data, we can drop all goal-keepers and repeat PCA."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first create a new data-frame with all goal-keepers removed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df2 = df[df[\"Position\"] != \"GK\"].copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we remove all the columns associated with the attributes that are mostly associated with goal-keepers. We also remove the columns with `pc1` and `pc2`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df2 = df2.drop(['GK diving',\n",
    " 'GK handling',\n",
    " 'GK kicking',\n",
    " 'GK positioning',\n",
    " 'GK reflexes','pc1','pc2'],1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Repeat all the steps from the previous analysis: shift to zero-mean, obtain svd, plot explained variances."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# clear\n",
    "\n",
    "Y = df2.iloc[:,:-4].copy()\n",
    "\n",
    "B = Y - Y.mean()\n",
    "u, s, vt = np.linalg.svd(B, full_matrices = False)\n",
    "v = vt.T\n",
    "\n",
    "variance = s**2\n",
    "sum_var = sum(variance)\n",
    "var_exp = [vv/sum_var for vv in variance]\n",
    "cum_var = np.cumsum(var_exp)\n",
    "\n",
    "plt.figure(figsize=(10,6))\n",
    "plt.bar(range(29),var_exp,label='individual explained variance')\n",
    "plt.plot(range(29),cum_var,'ro-',label='cumulative explained variance')\n",
    "plt.legend(loc=0)\n",
    "plt.xlabel(\"components\")\n",
    "plt.ylabel(\"% variance\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add the first two components to the data frame and plot them in a scatter plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# clear\n",
    "df2['pc1'] = u[:,0]*s[0]\n",
    "df2['pc2'] = u[:,1]*s[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#clear\n",
    "\n",
    "g = sns.lmplot(x = \"pc1\", y = \"pc2\", data = df2, hue = \"Position\", fit_reg=False, height=11, aspect=2, legend=True,\n",
    "            markers=[\"o\", \"^\",\"s\"],palette=dict(FWD=\"g\",  MID=\"r\", DEF=\"m\"))\n",
    "ax = g.axes[0,0]\n",
    "ax.axvline(x=0,color='k', ls = '--')\n",
    "ax.axhline(y=0,color='k', ls = '--')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the weights of each attribute corresponding to the principal component 1:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#clear\n",
    "\n",
    "labels_new = Y.columns.values\n",
    "\n",
    "plt.figure(figsize=(14,6))\n",
    "plt.bar(labels_new,v[:,0])\n",
    "plt.xticks(rotation=90);\n",
    "plt.title('importance of each attribute in ${\\\\bf p}_1$');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " In the same figure, plot the projections for the attributes in columns [4,9,11,12,14,3,12] onto principal components 1 and 2. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#clear\n",
    "\n",
    "g = sns.lmplot(x = \"pc1\", y = \"pc2\", data = df2, hue = \"Position\", fit_reg=False, height=11, aspect=2, legend=True,\n",
    "            markers=[\"o\", \"^\",\"s\"],palette=dict(FWD=\"g\",  MID=\"r\", DEF=\"m\"))\n",
    "ax = g.axes[0,0]\n",
    "ax.axvline(x=0,color='k', ls = '--')\n",
    "ax.axhline(y=0,color='k', ls = '--')\n",
    "\n",
    "s = 300   # this will scale the size of the arrow plot\n",
    "\n",
    "for J in [4,9,11,12,14,3,12]:\n",
    "\n",
    "    x = v[J,0]    \n",
    "    y = v[J,1]  \n",
    "\n",
    "    # make an arrow from the origin to a point at (x,y)\n",
    "    ax.arrow(0,0,s*x,s*y,color='black',width=1)\n",
    "    ax.text(x*s*1.5,y*s*1.1,labels_new[J],fontsize=28)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#clear\n",
    "\n",
    "g = sns.lmplot(x = \"pc1\", y = \"pc2\", data = df2[df2['Position'] == 'DEF'], hue = \"Position\", fit_reg=False, height=11, aspect=2, legend=True,\n",
    "            markers=[\"s\"],palette=dict(DEF=\"m\"))\n",
    "ax = g.axes[0,0]\n",
    "ax.axvline(x=0,color='k', ls = '--')\n",
    "ax.axhline(y=0,color='k', ls = '--')\n",
    "\n",
    "\n",
    "s = 300   # this will scale the size of the arrow plot\n",
    "\n",
    "for J in [4,9,11,12,14,3,12]:\n",
    "\n",
    "    x = v[J,0]    \n",
    "    y = v[J,1]  \n",
    "\n",
    "    # make an arrow from the origin to a point at (x,y)\n",
    "    ax.arrow(0,0,s*x,s*y,color='black',width=1)\n",
    "    ax.text(x*s*1.5,y*s*1.1,labels_new[J],fontsize=28)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 3: Food consumption in UK (using sklearn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Example from:\n",
    "http://setosa.io/ev/principal-component-analysis/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.io.parsers.read_csv(\"additional_files/UK_foods.csv\")\n",
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/Picture1.png\" width=\"500\"/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "headers = data['Unnamed: 0'].values.tolist()\n",
    "print(headers)\n",
    "\n",
    "new_data = data.drop(['Unnamed: 0'], axis=1)\n",
    "new_data.head()\n",
    "\n",
    "regions = new_data.columns.values.tolist()\n",
    "print(regions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "food = pd.DataFrame(new_data.values.T,columns=headers)\n",
    "food['region'] = regions\n",
    "food"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "food.std()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### This is when we want to try PCA!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Performing PCA without scaling the data (to match the results from the website)\n",
    "#X = pd.DataFrame(food[headers], columns=headers)\n",
    "# In general, PCA scales the variables to zero-mean (use line below to scale)\n",
    "X = pd.DataFrame(scale(food[headers]), columns=headers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca = PCA().fit(X)\n",
    "pca_samples = pca.transform(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var_exp = pca.explained_variance_ratio_\n",
    "plt.bar(range(len(var_exp)),var_exp, align='center', label='individual explained variance');\n",
    "plt.ylabel('Explained variance ratio');\n",
    "plt.xlabel('Principal components');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "components = pd.DataFrame(pca.components_, columns = headers) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.bar(headers,components.values[0])\n",
    "plt.xticks(rotation=90)\n",
    "plt.title('influence of original variables(food) upon pc1')\n",
    "plt.figure()\n",
    "plt.bar(headers,components.values[1])\n",
    "plt.xticks(rotation=90)\n",
    "plt.title('influence of original variables(food) upon pc2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xstar = pd.DataFrame(pca_samples,columns=['pc1','pc2','pc3','pc4'])\n",
    "Xstar['region'] = regions\n",
    "Xstar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.stripplot(x=\"pc1\",y=\"region\", data=Xstar, jitter=0.05, linewidth=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = plt.figure()\n",
    "ax = sns.lmplot('pc1', 'pc2',Xstar,hue='region', fit_reg=False)\n",
    "plt.axis('equal')\n",
    "plt.xlabel('pc1')\n",
    "plt.ylabel('pc2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
